Black hole-neutron star mergers: effects of the orientation of the black hole spin 



(N, 

^: 
^: 

^': 

Oh- 
I ' 

^ . 



> 

m 

o 

^, 

1^ 

O 

o 



X 



Francois Foucart,' Matthew D. Duez,'-^ Lawrence E. Kidder,' and Saul A. Teukolsky''-' 

' Center for Radiophysics and Space Research, Cornell University, Ithaca, New York, 14853, USA 
^ Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA 
Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, California 91125, USA 

The spin of black holes in black hole-neutron star (BHNS) binaries can have a strong influence on the merger 
dynamics and the postmerger state; a wide variety of spin magnitudes and orientations are expected to occur in 
nature. In this paper, we report the first simulations in full general relativity of BHNS mergers with misaligned 
black hole spin. We vary the spin magnitude from aBn/A/BH = to aBn/A/BH ~ 0.9 for aligned cases, and 
we vary the misalignment angle from to 80° for obh/A/bh ~ 0.5. We restrict our study to 3:1 mass ratio 
systems and use a simple F-law equation of state. We find that the misalignment angle has a strong effect on the 
mass of the postmerger accretion disk, but only for angles greater than ^ 40° . Although the disk mass varies 
significantly with spin magnitude and misalignment angle, we find that all disks have very similar lifetimes 
~ 100ms. Their thermal and rotational profiles are also very similar. For a misaligned merger, the disk is tilted 
with respect to the final black hole's spin axis. This will cause the disk to precess, but on a time scale longer 
than the accretion time. In all cases, we find promising setups for gamma-ray burst production: the disks are 
hot, thick, and hyperaccreting, and a baryon-clear region exists above the black hole. 

PACS numbers: 04.25.dg, 04.40.Dg, 04.30.-w, 47.75.-l-f, 95.30.Sf 



I. INTRODUCTION 

Black hole-neutron star (BHNS) binary mergers present a 
remarkable opportunity to study strongly-curved spacetime 
and supernuclear-density matter in the most extreme, dynam- 
ical conditions. BHNS binaries in compact orbits emit strong 
gravitational waves, and they are expected to be one of the 
main sources for Advanced LIGO and VIRGO ||il|2|l- Current 
estimates for the event rates of binary mergers coming from 
population synthesis models predict that Advanced LIGO will 
see about 10 BHNS/yr, although uncertainties in the models 
allow for a large range of potential event rates, ^ 0.2 - 300 
BHNS/yr |3]. These gravitational waves contain, in principle, 
a wealth of information on their source, such as the mass and 
spin of the black hole (BH) and the mass and radius of the 
neutron star (NS). Inferred properties of the NS could be used 
to constrain the NS equation of state. Information in the waves 
can only be extracted, however, by comparison with accurate 
numerically-generated predictions that provide the expected 
waveform for each possible BHNS system. 

Mergers of BHNS binaries have also been proposed 
as potential progenitors of short-hard gamma-ray bursts 
(SGRB) lH. The origin of SGRBs is not yet known, although 
it is certain that the engines are compact and located at cosmo- 
logical distances, and there is evidence (such as their presence 
in nonstar forming regions) to support a mechanism differ- 
ent from that associated with long-soft GRBs, namely stellar 
core collapse. For BHNS mergers, the generation of a SGRB 
is possible only if the remnant black hole is surrounded by a 
massive, hot, thick accretion disk. Also, to obtain relativis- 
tic jets and a beamed outflow, a region mostly devoid of any 
matter is necessary (see e.g. jsll and references therein). Only 
numerical simulations in full general relativity with realistic 
microphysics can determine if these conditions are likely to 
be obtained. 

Whether a disk forms or not will depend on the premerger 
characteristics of the binary, especially the BH mass, the NS 



radius, and the BH spin. Current estimates from population 
synthesis models suggest that most systems are likely to be 
formed with a black hole of ^ 10 A/©. Relativistic simulations 
to date have considered cases of relatively low mass black 
holes (~ 2 - 7Mq ) 16^1311, for which the NS is expected 
to disrupt outside the innermost stable circular orbit (ISCO), 
making disk formation more likely. These simulations have 
found cases of massive disk formation, with A/bh ~' 3 — AMq 
resulting in the largest disks II12I1 . The NS radius is the param- 
eter related to the equation of state that has the largest effect 
on the waveform and post- merger disk lll3ll . with larger radii 
resulting in larger disks JSl llSll . 

The spin of the black hole can have a strong influence on 
the merger. The ISCO is smaller for prograde orbits around a 
spinning BH than for orbits around a nonspinning hole. Be- 
cause disk formation is expected to be more likely if NS tidal 
disruption occurs outside the ISCO than if it occurs inside, 
BH spin can facilitate disk formation. With high BH spin, 
it is even plausible that BHNS binaries with the most likely 
mass ratios (~7:1) give rise to substantial disks UM- The 
magnitude of the BH spin is largely unconstrained by pop- 
ulation synthesis models, as it comes mostly from the spin 
acquired during formation of the black hole in a core-collapse 
event II15I1 . The effect of aligned and antialigned spins was in- 
vestigated in full general relativity for the 3:1 mass-ratio case 
by Etienne et al. u2i\ . They found that a large aligned spin 
(and correspondingly small ISCO) leads to a much more mas- 
sive post-merger disk. For example, for aBnA'^BH = 0.75 
and A/bh ~ 4Afo, a disk of A/disk ^ O.2Af0 can be ob- 
tained. For BHNS binaries with massive black holes (A/bh ^ 
10A/q), forming a disk may in fact only be possible if the hole 
is spinning. 

There is no reason to expect the black hole spin to be 
aligned with the orbital angular momentum. Population syn- 
thesis models predict a relatively wide distribution of orien- 
tations, with about half of the binaries having a misalign- 
ment between the BH spin and the orbital angular momen- 



turn of less than 45° when the initial BH spin is abh A^bh = 
0.5 iflSll . Misalignment can reduce or reverse the BH spin 
effects described above. This can be understood by consider- 
ing prograde orbits of test particles with small radial velocity, 
which become unstable farther away from the BH for inclined 
orbits than for equatorial orbits of the same angular momen- 
tum. Misalignment will also produce qualitatively new ef- 
fects, including the precession of the premerger orbital plane 
and BH spin. The influence of misalignment has been stud- 
ied for 10:1 mass-ratio binaries in the approximation that the 
spacetime is assumed to be Kerr. Rantsiou et al. Illql showed 
that, in this approximation, disks can be formed only at rela- 
tively low inclinations and only for near extremal black holes. 
Ultimately, though, simulations in full general relativity are 
needed to accurately model such systems. 

In this paper, we report on fully relativistic studies of mis- 
aligned BHNS binaries. We limit ourselves to small mass sys- 
tems (^/bh ~ 4.2Mq) and a simplified equation of state, 
but we consider a significant range of black hole spin mag- 
nitudes and orientations. We confirm the results of Etienne et 
al urn regarding the effects of an aligned BH spin. For mis- 
aligned spins, we find that the misalignment angle can have 
a strong effect on the post-merger disk mass, but only for an- 
gles greater than around 40°. Although the disk mass varies 
greatly with BH spin, most other disk properties are very simi- 
lar, including the accretion time scale, the location of the max- 
imum density, the average temperature, and the entropy and 
angular momentum profiles. The disks are all thick, each with 
a height H to radius r ratio H/r sa 0.2, nearly independent 
of r. Crucially, they all have a baryon-clear region above and 
below the BH. The disks are misaligned with the final BH spin 
by < 15°. They do precess about the BH spin axis, but with- 
out reaching a fixed precession rate. Indeed, the steady-state 
precession time scale is expected to be significantly longer 
than the accretion time scale. 

This paper is organized as follows. In Sec. [Ill we discuss the 
method used to construct the very general BHNS initial data 
we use. We also discuss in detail the improvements to our 
evolution code that have increased the accuracy by an order 
of magnitude over the results presented in Duez et al. um . 
We then present our run diagnostics in Sec. |III] The different 
cases to be evolved are described in Sec.|IVl We then present 
the results of the simulations in Sees. |V] and [Vl] Finally, we 
draw conclusions in Sec. I VIII 



n. NUMERICAL METHODS 

A. Initial data 

For numerical evolutions of Einstein's equations, we de- 
compose the spacetime under study into a foliation of space- 
like hypersurfaces parametrized by the coordinate t. Ein- 
stein's equations can be written in the form of hyperbolic evo- 
lution equations plus a set of constraints that have to be satis- 
fied on each t ~ constant slice. Our initial data at t = must 
be chosen such that it satisfies these constraints. We construct 
initial data using the Extended Conformal Thin Sandwich for- 



malism (XCTS) lUTl \\m . If we write the spacetime metric 
as 

ds — g^^dx^dx'^ 

= ^a^dt^ + V^7.y(dx* + I3'dt)idx^ + j3^dt), (1) 

the initial data to be determined include the lapse a, the shift 
vector /?', the conformal factor ip, the conformal 3-metric 7^^ 
and the extrinsic curvature /v,,^ = —\Cngi_Lu (where £„ is 
the Lie derivative along the normal n to the f = slice). 
The constraints can be expressed as a set of 5 coupled ellip- 
tic equations for the lapse, shift, and conformal factor lUSll . 
The physical properties of the system are then determined by 
the choice of the remaining free parameters: the trace of the 
extrinsic curvature K = g^^^K^^'^, the conformal metric ■jij, 
their time derivatives dtK and dtjij, and the matter stress- 
energy tensor T'^]^*-^^'^. 

The system of elliptic equations is solved using the spec- 
tral elliptic solver SPELLS developed by the Cornell-Caltech 
collaboration Ill9ll . and initially used to construct initial data 
for binary black hole systems by Pfeiffer et al. lllSl uOT . A 
detailed presentation of the methods used for the construction 
of BHNS initial data was given in Foucart et al. 112111 . Here, 
we limit ourselves to a brief summary plus a description of 
the changes made to accommodate the possibiUty of arbitrary 
spin orientation. 

As the system is expected to be initially in a quasiequilib- 
rium state, with the binary in a low-eccentricity circular orbit 
of slowly decreasing radius, we work in a frame comoving 
with the binary and set the time derivatives to zero: dtK = 
and dtjij = 0. As for 7^ and K, we make a choice inspired 
by the results of Lovelace et al. 1I22I1 for binary black hole 
systems. Close to the BH, the metric matches its Kerr-Schild 
values for a BH of the desired mass and spin, while away from 
the BH, the conformal metric is flat and K = 0. The transi- 
tion between these two regions is done by using the following 
prescription: 



K 
X 

VBH 



^ij + [lii (^BH, vbh) - Sij]e' 



->^{r-Bn/w) 



= K *(aBH,VBH)e 

J'BH - ?'AH 



-Kr-Ba/w) 



0'°* X CBH, 



, (2) 
(3) 
(4) 
(5) 



where the KS subscript refers to the Kerr-Schild values, 
''BH (''Ns) is the coordinate distance to the center of the BH 
(NS), rAH the coordinate radius of the BH apparent hori- 
zon, obrZ-^bh is the dimensionless spin of the BH, q ^ 
A/ns/^'^bh a constant of the order of the mass ratio, cbh the 
coordinate location of the BH center with respect to the center 
of mass of the system, and w is some freely specifiable width, 
chosen so that the metric is nearly flat at the location of the 
NS. The parameter A, which is designed to impose a flat back- 
ground at the location of the NS, is set to 00 for tns < 9?'bh- 
Boundary conditions are imposed at infinity and on the ap- 
parent horizon of the BH (since the inside of the BH is excised 
from our computational domain). The boundary conditions at 
infinity are chosen so that the metric is asymptotically flat. 



while the inner boundary conditions follow the prescriptions 
of Cook and Pfeiffer |[23ll . which make the inner boundary 
an apparent horizon in quasiequilibrium. There is some free- 
dom in these boundary conditions: on the apparent horizon, 
the conformal lapse is not fixed (we set it to the value of an 
isolated Kerr BH), and the shift is determined only up to a ro- 
tation term P'' = /3* + e'^'^nf^Xk- The value of fl^^ deter- 
mines the spin of the BH, but the exact relation between H^^ 
and the spin is unknown a priori; to get the desired BH spin, 
we have to solve iteratively for Jl^^. On the outer boundary, 
the shift can be written as 



/3 = n™' X r + aor + t-boost 



(6) 



where ff^"^ allows for a global rotation of the coordinates, a'o 
for a radial infall with velocity v = aor, and «boost for a 
boost. As an initial guess for the orbit of the binary, we can 
set the radial velocity at t = to (a'o = 0). This assumption, 
as well as the quasiequilibirum formalism, clearly neglects the 
evolution of the orbit over time through the radial infall of the 
binary and the precession of the orbital plane. Both effects 
are, however, acting over relatively long time scales: over its 
first orbit, even the binary with the most inclined spin consid- 
ered here (s.5i80 in the later sections) goes through less than 
10% of a full precession period of the BH spin while the coor- 
dinate separation between the compact objects is reduced by 
about 20%. One known effect of the quasi-circular approxi- 
mation is that the binary will have a nonzero eccentricity. The 
eccentricity can be decreased bymodifying the initial values 
of a'o and $1'°' f24fl (see also |21|| for an application of that 
method to BHNS binaries) as long as the initial eccentricity 
and orbital phase can be accurately measured. Here, we only 
apply this technique when the spin of the BH is aligned with 
the orbital angular momentum of the binary. For precessing 
binaries, significantly reducing the eccentricity would require 
a larger initial separation for which the effects of eccentricity, 
precession and radial infall can be properly disantangled. 

In the presence of matter, additional choices are required. 
We assume that the fluid is in hydrostatic equilibrium in the 
comoving frame, and require that its 3-velocity is irrotational. 
The first condition gives an algebraic relation between the en- 
thalpy h of the fluid, its 3-velocity Vi, and the metric g^^, 
while the second leads to another elliptic equation determin- 
ing the velocity field. These equations are coupled to the con- 
straints: the whole system can only be solved through an it- 
erative method. For a BH with a spin aligned with the total 
angular momentum of the system, that method is described in 
Foucart et al. UM : we solve for the metric using SPELLS, then 
determine the new value of the enthalpy h, as well as the or- 
bital angular velocity Jl™' (chosen so that the binary is in qua- 
sicircular orbit), the position of the BH in the equatorial plane 
(so that the total linear momentum Padm vanishes), and the 
free parameter l^gjj (to drive the spin of the BH to its desired 
value). Finally, we solve for the velocity field through the el- 
liptic equation imposing an irrotational configuration, and go 
back to the first step. 

In order to construct initial data for BHs with a spin that is 
not aligned with the orbital angular momentum of the binary, 
a few changes are necessary. First, we do not assume that 
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FIG. 1 : Evolution of the position of the center of mass along the 
z-axis during inspiral (solid line), compared to the motion expected 
from the boost v^ given in the initial data (dotted line). 



rj^^ is aligned with the orbital angular momentum. Instead, 
all 3 components of $7^^ are solved for. We also abandon the 
assumption of equatorial symmetry, and control the position 
of the BH along the z-axis of the orbital angular momentum 
by requiring that the NS is initially moving in the x^^-plane, 
with its center in the 2 = plane (the z coordinate of the lo- 
cation of the BH is chosen so that the condition e^-Vh ~ 
is satisfied at the center of the NS). Finally, to guarantee that 
^ADM = we impart a boost to the whole system through the 
boundary condition at infinity: v^^^^^ = v^. The center of 
mass is then expected to have a global motion during inspiral 
corresponding to that boost, and we check in Fig. [T] that this 
is indeed the case. By adding these conditions to the iterative 
procedure used to generate BHNS initial data, we are able to 
obtain high-precision initial configurations for arbitrary val- 
ues of the orientation of the BH spin. 



B. Evolution 

The simulations presented here use the SpEC code devel- 
oped by the Cornell-Caltech-CITA Collaboration iH. To 
evolve BHNS systems, the two-grid method described in Duez 
et al. lUir is used. Einstein's equations are evolved on a pseu- 
dospectral grid, using the first-order generalized harmonic for- 
mulation 12a], while the hydrodynamical equations are solved 
on a separate finite difference grid called the "fluid grid". The 
hydrodynamic equations are written in conservative form 



atU + vF(u) = s(u). 



(7) 



To compute the flux F on the faces of each finite difference 
cell, we use the third-order shock capturing PPM reconstruc- 
tion method 12711 . More details on the numerical methods used 
can be found in Duez et al. UAn . However, since the publica- 
tion of |ill], several important improvements have been made 
to the code, which are described in the following subsections. 



1. Dynamic regridding 



region far away from it. The exact form of the map is 



To accurately evolve a BHNS binary while determining 
its gravitational wave emission, simulations have to resolve 
events occurring at very different scales. When the neutron 
star is disrupted and a disk forms, we expect shocks in the 
disk, and steep density and temperature gradients close to the 
BH. But the disruption of the star also leads to the creation of 
a long tidal tail which can initially contain up to 5 — 10% of 
the initial mass of the star and expand hundreds of kilometers 
away from the center of the BH II13I1 . Clearly, both the sharp, 
small-scale features around the black hole and the large-scale 
tidal tail should be properly resolved if we want to follow the 
formation of an accretion disk. Furthermore, to extract grav- 
itational waves accurately, the evolution of the gravitational 
fields should extend to the wave zone, in regions where no 
matter at all is present. 

Because we use different grids to evolve the metric and the 
fluid variables, the spectral grid on which we solve the gen- 
eralized harmonic equations can be extended into the wave 
zone while the fluid grid used for the relativistic hydrodynam- 
ical equations only covers the region where matter is present. 
Our earlier simulations lUlll were limited to nonspinning black 
holes. In that case, most of the matter was rapidly accreted 
onto the hole and the tidal tails and accretion disks were small 
enough that manually expanding the fluid grid at a few cho- 
sen timesteps allowed us to resolve the evolution of the fluid 
at a reasonable computational cost. For spinning black holes, 
this is no longer the case: cost-efficient evolutions require a 
grid with points concentrated in the high-density regions close 
to the BH, and coarser resolution in the tail. Furthermore, 
as the evolution of the fluid is highly dynamical, interrupting 
the simulation whenever the finite difference grid is no longer 
adapted to the fluid configuration becomes impractical. One 
solution would be to use an adaptive mesh refinement scheme, 
similar to the codes used by Yamamoto et al. ||28I1 and Eti- 
enne et al. ||12ll. In our code, we choose instead to use a map 
between the fluid grid and the pseudospectral grid that con- 
centrates grid points in the region close to the black hole and 
automatically follows the evolution of the fluid. 

To do this, we measure the outflow of matter across surfaces 
close to but inside of the fluid grid boundaries. As soon as the 
outflow across one of these surfaces crosses a given threshold 
(chosen so that the amount of matter leaving the grid over the 
whole simulation is negligible compared to the final mass of 
the accretion disk), the grid expands. The opposite is done on 
fixed surfaces farther away from the grid boundaries to force 
the grid to contract whenever the fluid is moving away from a 
boundary. The map itself is the combination of: 

(i) A translation of the center of the grid, to follow the gen- 
eral motion of the fluid 

(ii) A linear scaling of each coordinate axis, to adapt to its 
expansions and contractions 

(iii) A radial map smoothly transitioning from a high resolu- 
tion region close to the black hole to a lower resolution 



r <rA 

f{r) - f{rA) + TA, rA <r < rs 

K^ -rs) + f{rB) - .f{rA)+rA, r > rs 

(8) 



/(r) = r{ar^ + br'^ +cr + d), 



(9) 



where the coefficients (a, b, c, d) are chosen so that the 
map is C^ at rA and rs and A is chosen so that the grid 
is of the desired size. The radii r^ and vb are fixed for 
the whole evolution and will determine respectively the 
minimum resolution in the neighborhood of the black 
hole and the characteristic lengthscale of the transition 
between the high and low resolution regions. 



2. Excision 

We find that our hydrodynamics code is more stable near 
the excision zone if we switch from PPM to the more dif- 
fusive MC reconstruction 02911 in the vicinity of the excision 
zone. Therefore, we replace the face values determined by 
PPM reconstruction, uh^l^^'^^ with a weighted average: 



UR,L ^ fUR,L'''^ + (1 - f)uR^ 



PPM 



(10) 



-[{r-ri)/rif f^^ 



where / = 1 for r < ri ^ 2rcx and / = 
r > ri. 

The MC face-value computation must be altered when its 
regular stencil would extend into the excised region, and do- 
ing this properly turns out to be important for stability. Con- 
sider a one-dimensional problem with grid points Xn = nAx. 
Then the face-value reconstruction of the function Ui from 
the left ML_i_i/2 = Wi-i/2-e and from the right Ufl,^j_i/2 = 
Ui-i/2+e must be adjusted as follows. 
(i) If Xi is in the excision zone, but Xi^i is outside, 
set UL,j+i/2 = "i-i and Ufl,i+i/2 = u^-i 
(ii) If Xi is in the excision zone, but Xi+i is outside, 
set UL,i-i/2 = Ui+i and Ufl,i_i/2 = Uj+i 
(iii) If Xi is outside the excision zone, but Xi^i is inside, 

SetUi_i_l/2 =MH^i-l/2 

(iv) If Xi is outside the excision zone, but Xi+i is inside, 

set Uji,i+i/2 = Ul,i+1/2 

We also observed that the stability of our code close to the 
excision surface was strongly affected by the details of the in- 
terpolation method chosen for the communication from finite 
difference to spectral grid in that region. Previously, the inter- 
polation stencil was shifted away from the excision boundary 
until the entire stencil was out of the excision zone. This could 
lead to unstable evolutions or large interpolation errors if the 
excision region happened to be located close to the boundary 
between two subdomains of the finite difference grid, and ac- 
ceptable stencils could only be found far from the point we 
were interpolating to — or could not be found at all (to limit 
MPI communications, the stencil has to be entirely contained 
in one subdomain). Currently, we limit the displacement of 



the stencil to a maximum of 3 grid point separations. If there 
is no good stencil within that distance, we decrease the or- 
der of the interpolation, and keep doing so until an acceptable 
stencil is found. 

Another interpolation method would be to forbid any dis- 
placement of the stencil, and immediately drop to lower order 
as soon as part of the stencil lies within the excision zone. 
Both algorithms are equally robust, but when tested on an ac- 
tual BHNS merger the first appeared to perform better at main- 
taining a smooth solution and low constraint violations on the 
excision surface. Accordingly, we chose it as our standard in- 
terpolation method and used it for all simulations presented in 
this paper. 
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3. Coordinate evolution 

In the generalized harmonic formulation, the evolution of 
the inertial coordinates Xa is given by the inhomogeneous 
wave equation 



W'WbXa = Ha 



(11) 



where Vf, is the covariant derivative along Xb- The evolution 
of the function Ha{xb) can be freely specified, but its value 
on the initial slice t = is determined by the initial data (the 
lapse and shift at t = fix the initial evolution of the gauge). 
While the binary spirals in, we choose dtHa{t, i,) =0 in the 
coordinate frame Xi comoving with the system. In our pre- 
vious paper jllll . we changed the gauge evolution during the 
merger phase by damping Ha exponentially in the comoving 
frame: 



^a v^7 *^i 



-(*-td)/T 



^a {J-d-) -^i 



(12) 



where td is the disruption time — the time at which we begin 
damping — and r is a damping time scale of order 1071/ {M 
being the total mass of the system). Further experimentation 
has shown that it is better not to change Ha near the excision 
zone. In our current simulations, we set 



Ha{t,Xi) 
Ha{td,Xi) 

Q{r) 



{Q(f) + [l-Q(f)]e-(*-*'')/^} (13) 

(14) 



^ g-('=/?cx)' + l 



during the merger phase, where r is the distance to the cen- 
ter of the black hole in the comoving frame, and fex is the 
excision radius. 



FIG. 2: Average surface density of the fluid at a given distance from 
the center of mass of the system, for the simulation s.5i20 described 
in Sec. lIVI The surface density is plotted at time io = fmorgcr + 
10.3ms, when we begin to evolve the system using the fixed-metric 
approximation, as well as 0.5ms and I.Ims later. We see that the 
profile is very similar for both evolution methods, even though the 
disk itself is not in a stationary configuration. 



a strong influence on the behavior of the system. In numeri- 
cal simulations, we can thus extract some information on the 
long-term behavior of the final black hole-accretion disk sys- 
tem by neglecting the evolution of the metric and only evolv- 
ing the fluid (c.f. 13011 ). Using this approximation, our code 
runs about 4 times faster 

To test the limitations of this method, we evolve the cou- 
pled system ^ 1ms past the time at which we begin the ap- 
proximate, fluid-only evolution. We look for differences in the 
accretion rate or the characteristics of the disk (temperature, 
density, inclination) between the two methods used. As long 
as we wait for the properties of the black hole to settle down 
before switching to the approximate evolution scheme, the 
two methods show extremely good agreement — except for 
the highest spin configuration, which leads to a massive disk 
that cannot be evolved accurately by fixing the background 
metric. In Fig. |2] we show the evolution of the density profile 
of a disk using both the full GR evolution and the fixed-metric 
approximation. The evolution of the accretion disk is mostly 
unaffected by the change of evolution scheme. 



III. DIAGNOSTICS 



4. Evolutions with fixed metric 

During the merger of a BHNS binary, both the spacetime 
metric and the fluid configuration are highly dynamical. Ein- 
stein's equations have to be solved together with the conserva- 
tive hydrodynamics equations, and the evolution of that cou- 
pled system is computationally intensive. However, a few mil- 
liseconds after merger, the BH remnant settles into a quasis- 
tationary state as it accretes slowly from the surrounding ac- 
cretion disk. Then, the evolution of the metric does not have 



An indispensable test of numerical accuracy is convergence 
with grid resolution. We have evolved most of the cases dis- 
cussed below at three resolutions. We call these Resl, Res2, 
and Res3; they correspond to 100"^, 120"^, and 140'^ gridpoints, 
respectively, on the fluid grid and to 69'^, 82'^, and 95'^ collo- 
cation points on the pseudospectral grid. 

The black hole is described by its irreducible mass 
A/jir, its spin Sbh, and its Christodoulou mass 7\/bh = 
^M^j, + S'bh^/(4M^j,). The spin Srh is computed using 
the approximate Killing vector method ll22ll . 



To monitor the nuclear matter, we first measure the bary- 
onic mass i\/b on the grid as a function of time. Initially, this 
will be the baryonic mass of the neutron star After the tidal 
destruction of the star, it will be the sum of the accretion disk 
and the tidal tail masses. At late times (more than 10ms after 
merger), it will be the baryonic mass of the disk. The accre- 
tion time scale Tdisk is Mb/{dMh/dt). 

We also analyze the heating in the disk. This is done 
through both an entropy and a temperature variable. When 
evolving with microphysical equations of state, the tempera- 
ture and entropy are provided directly. For this study, we use 
a F-law equation of state, so we need a way to estimate the 
entropy and temperature. The entropy measure s we define as 
s = \og{Hi/ K,i), where k is the polytropic constant obtained 
from the relation between the pressure and the baryon density 
{P = Kpo)' ^^^ ^i is i'^s initial, cold value. To estimate the 
physical temperature, T, we assume that the thermal contri- 
bution to the specific internal energy eth = e ^ ^{T = 0) is 
given by a sum of ideal gas and blackbody components: 



angular momentum Jdisk are given by 



Eth 



3fcT 
2m„ 



,aT'^ 



(15) 



where m„ is the nucleon mass, and the factor / reflects the 
number of relativistic particles, and is itself a function of T. 
(See ||all2l. who also make this assumption.) For the most 
part, we will report density-averaged values of s and T. For 
example, the density-averaged entropy is 



L^"' = 



S" 



(.t'^T'^o - xT^°) d^x 
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These parameters determine the inclination and the preces- 
sion of the disk with respect to the spin of the black hole: if 
Jrh ~ JbrGz, then the orbital angular momentum of the disk 
at radius r is written as 

Jdisk (j') = >/disk (t") (sin /? cos 76^; + sin /3 sin jCy + cos iSe^) . 

(22) 
Another useful property of the disk is its scale height, H. 
For a disk with exponentially decreasing density, H is defined 
by p ~ pce~^/^ . Here, however, the vertical profile of the 
disk is significantly more complex, and various definitions of 
H could be considered. We use the spread of the density dis- 
tribution p{9, (p) on a sphere of constant radius r with its polar 
axis along Jdisk (j') and define 



H(r) = r tan 



Jpir)dS _ 



(23) 
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To launch a GRB, a baryon-clean region above the disk is 
probably needed. This does not mean that a wider clean region 
is always better, since a thick disk can help collimate the out- 
flowing jet — but we want to determine whether such a region 
exists or not. To estimate the baryon-poor opening above our 
disks, we define the opening angle 6'cican- This angle specifies 
the widest cone oriented along the BH spin in which the con- 
dition p < pcut is everywhere satisfied: if ^cioan('', 4') is the 
opening angle within which we have p < pcut at radius r and 
azimuthal coordinate cj), then ^cican ~ minr.0 [6'cican (?', 0)]- In 
these simulations, the numerical method requires atmospheric 
corrections to be applied starting at p = 6 x lO^g/cm'^, so 
that it is impossible to reliably predict the behavior of matter 
below that threshold. Therefore, we set pcut = 3 x lO^g/cm'^. 

For precessing binaries we also compute the tilt B and twist 
7 of the disk, as defined by Fragile & Anninos II31I1 . If x'^ are 
the inertial coordinates, T'^'^ the stress-energy tensor, e^^cri 
the Levi-Cevita tensor (with its last index limited to nonzero 
values), Jbh the angular momentum of the BH and iy an ar- 
bitrary unit vector orthogonal to Jbh, then /3, 7 and the disk 



The parameter /i is somewhat arbitrary. For an exponential 
profile 1^1 ~ 0.5, while for a constant density profile [p = pa 
for 6 < H/r and p = otherwise) we have p ^ 3. The disks 
observed in our simulations are somewhat in between these 
two extremes. Accordingly, we make the approximate choice 

To measure the accuracy of our simulations, we monitor the 
ADM Hamiltonian and momentum constraints, and the gen- 
eralized harmonic constraints ||C|| i26ll . At our middle res- 
olution, ||C|| peaks below 1% for all cases and is less than 
0.1% during most of the inspiral. We also monitor the ADM 
mass il/ADM and angular momentum Jadm- An important 
check of our simulations is that the changes in these quanti- 
ties match the flux of energy and angular momentum in the 
outgoing gravitational radiation, which we reconstruct from 
the Newman-Penrose scalar V'4 as in Ref. 0211 . 



IV. CASES 

In order to assess the influence of the black hole spin on the 
disruption and merger of BHNS binaries, we study configura- 
tions for which all other physical parameters are held constant. 
The mass of the black hole is TI/bh = SMns, where 7\/ns 
is the ADM mass of an isolated neutron star with the same 
baryon mass as the star under consideration, and the initial 



coordinate separation is d ~ 7.5M, with M = Mbh + Mns- 
For the nuclear equation of state, we use the polytrope 



p = (r- i)p€ = Kp^ + Tp 



(24) 



where T is a fluid variable related to, but not equal to, the 
physical temperature. We set F = 2 and choose k so that 
the compaction of the star is C = M-^^s/Rns ~ 0.144. For 
poly tropic equations of state, the total mass of the system does 
not have to be fixed: results can easily be rescaled by M (see 
e.g. Sec. II-F of Foucart ef a/. IuHI ). However, whenever we 
choose to interpret our results in physical units (ms, km, Mq), 
we will assume that A/ns = IAMq (M = 5.6A/0). For that 
choice, the neutron star has a radius i?NS = 14.6km, and the 
initial separation is d = 63km. 

The different initial configurations and black hole spins 
studied are summarized in Table U We consider 3 dif- 
ferent magnitudes of the dimensionless spin aBnA^BH = 
(0, 0.5, 0.9), all aligned with the orbital angular momentum. 
Then, we vary the inclination angle (/)bh between the spin 
of the black hole and the initial angular velocity of the sys- 
tem, $1'°'. Considering that most BHNS binary systems 
are expected to have 0bh < 90° (Belczynski et al. lUSll ). 
with about half of the binaries at (/jbh < 40°, we choose 
(/>BH = (20°, 40°, 60°, 80°). The orientation of the com- 
ponent of the BH spin lying in the orbital plane could also 
have measurable consequences. For example, Campanelli et 
al. 1I33II showed that the superkick configuration found in bi- 
nary black hole systems is sensitive to the direction of the mis- 
aligned component of the BH spin. For BHNS binaries, kicks 
are relatively small, and we are more interested in the charac- 
teristics of the final black hole-disk system. After looking at 
different orientations for (/)bh = 80°, we find that the influ- 
ence of the orientation of the misaligned component of the BH 
spin is negligible compared to the influence of 0bh- For this 
first study of misaligned spins, we will thus limit ourselves to 
configurations for which the initial spin lies in the plane gen- 
erated by the initial orbital angular momentum and the line 
connecting the two compact objects. As the different initial 
configurations do not use the same background metric, there 
is no guarantee that two binaries with the same initial coordi- 
nate separation can be directly compared. A better compari- 
son between initial configurations is the orbital angular veloc- 
ity of the system. In Table U we show that all configurations 
have initial angular velocity within 1% of each other. This 
is the level of error that we expect from the quasiequilibrium 
method for binaries at this separation. 



V. THE NONSPINNING CASE: A TEST OF OUR 
ACCURACY 

As an example, we consider the case sO, in which the BH is 
initially nonspinning. We evolve this case at each of our three 
resolutions. After a short (two orbits) inspiral, the neutron star 
is disrupted, and most of the matter is quickly swallowed by 
the black hole. The remainder expands into a tidal tail and 
then falls back to form an accretion disk. 



Case 


asn/AfBH 


0BH 


^initM 


-Efc/A/ADM 


Jadm/Maum 


C merger 


sO 





- 


4.16e-2 


9.5e-3 


0.66 


7.5ms 


s.5iO 


0.5 





4.11e-2 


l.Ole-2 


0.91 


11.4ms 


s.9i0 


0.9 





4.13e-2 


9.6e-3 


1.13 


15.0ms 


s.5i20 


0.5 


20 


4.09e-2 


l.Ole-2 


0.90 


10.5ms 


s.5i40 


0.5 


40 


4.10e-2 


l.Ole-2 


0.87 


9.9ms 


s.5i60 


0.5 


60 


4.11e-2 


9.9e-3 


0.82 


9.0ms 


s.5i80 


0.5 


80 


4.13e-2 


9.6e-3 


0.76 


7.7ms 



TABLE I: Description of the cases evolved. aBH/A/sH is the ini- 
tial dimensionless spin of the BH, ^bh is its inclination with respect 
to the initial orbital angular momentum and Eb is the initial binding 
energy, imcrgcr is defined as the time by which half of the matter has 
been accreted by the BH. Differences in the initial angular velocity 
and binding energy are within the margin of error of the initial data: 
at this separation the eccentricity reduction method can require vari- 
ations of fiinit of order 1%, and modifies the binding energy by a 
few percent. 



Nearly identical systems have been studied both by Shi- 
bata et al 131 and by Etienne et al 11211 . The former found 
an insignificant disk after merger, while the latter found 4% 
of the NS mass still outside the hole 300M (^ 8ms) after 
merger Both groups found a final BH spin of s = 0.56. We 
find a disk mass of 3.7% at 300A/ after merger, smaller than 
in 1 12], but closer to this result than to that in Igfl. Our final 
BH spin is 0.56, in agreement with both previous studies. 

In Fig. |3] we show the evolution of Mb and (s) for the en- 
tire merger phase for the three resolutions. Reassuringly, the 
different resolutions give very similar results, with the two 
higher resolutions being closest together. The baryonic mass 
is initially constant before accretion starts. Then, as the NS 
is disrupted and the core of the star is swallowed, Mb drops 
rapidly. It next levels off while the remaining matter is in an 
accretion disk and an expanding tidal tail. When the tidal tail 
falls back onto the disk, there is a second phase of rapid ac- 
cretion, after which the accretion rate settles down to a low 
value. At the end of the simulation, the accretion time scale 
is Tdisk ^ 55ms, implying that the total lifetime of the disk 
would be around 75ms. At late times, the deviation in Mb 
between resolutions becomes somewhat larger, indicating that 
our errors have accumulated to about 0.1% of the initial mass. 
For the purposes of this paper, this is adequate, since the disk 
on these time scales is affected by magnetic and radiation pro- 
cesses not included in the simulations. However, future long- 
term disk simulations will require higher accuracy. 

As for the entropy, at the beginning of the merger it only 
deviates from zero because of numerical heating during the 
inspiral. As expected, this numerical heating is significantly 
lower at higher resolutions. The post-merger heating is not nu- 
merical, but a physical consequence of shocks in the disk and 
the disk-tail interface. A confirmation that the heating is phys- 
ical is that it is nearly the same for all resolutions, and its mag- 
nitude is much larger than the numerical heating. When the 
disk settles, there is no further shock heating, so the entropy 
levels off. This indicates that the heating due to numerical 
viscosity is small compared to shock heating. Unfortunately, 
this is not the same as saying that the numerical viscosity is 




FIG. 3: Baryonic mass A/t normalized by its initial value Mi,,o and 
average entropy (s) for three resolutions. 
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FIG. 4: Madm and Jadm (normalized to their initial values) com- 
pared to the changes expected from the gravitational radiation flux 
for three resolutions. 



irrelevant altogether. However, the closeness of Tdisk at each 
resolution indicates that this viscosity is not the main driving 
force of the accretion. The average temperature (T) behaves 
in a way similar to the entropy. Starting from low values, it 
increases after the merger and stabilizes around 3MeV. All 
resolutions show the same (T) growth, and all level off at the 
same value. After leveling off, though, (T) displays O.lMeV 
oscillations that do not converge well, another indication that 
our accuracy is sufficient for some but not all purposes. 

In Fig. m we plot the ADM energy and orbital-axis angu- 
lar momentum measured on a surface 75i\/ from the center 
of mass of the system. Also plotted is the evolution of these 
quantities expected from the gravitational radiation through 
this surface. Overall, the agreement is quite good, although 
there is some deviation in A/adm a while after the merger. 
This seems to be associated with an increase in constraint 
violations at the merger time. The relative constraint viola- 
tions, as measured by j|C||, peak slightly below 1% at the mid- 
dle resolution. The corresponding values for the ADM con- 
straints are 1 -2 %, before both constraints fall back to low 
values. The deviations in Madm happen around the time this 
constraint-violating pulse reaches the r ~ 75M surface. 

We have checked several other quantities, including the 
black hole mass and spin and the gravitational waveform. All 
of these show very good convergence. 



VI. RESULTS 

The general behavior of our simulations is typical of BHNS 
binaries for which the NS is disrupted outside the innermost 
stable circular orbit of the BH. From an initial separation of 
60kni, the compact objects go through 2-3 orbits of inspiral 
driven by the emission of gravitational waves. When the dis- 
tance between them has been reduced to about 30 — 40km, 
tidal forces cause the neutron star to disrupt. Most of the mat- 
ter is rapidly accreted into the black hole, while the rest is di- 
vided between a long tidal tail, expanding about 200km away 
from the black hole, and a developing accretion disk. The 



Case 


Mdisk/A'/NS 


(T>disk 


/3(E^,x) 


r(E^ax) 


"clean 


— (^max) 


sO 


5.2% 


3.0MeV 


0° 


50km 


50° 


0.20 


s.5iO 


15.5% 


3.5MeV 


0° 


50km 


35° 


0.25 


s.9i0 


38.9% 


5.6 MeV 


0° 


20km 


8° 


0.18 


s.5i20 


14.5% 


3.6MeV 


2° 


50km 


40° 


0.20 


s.5i40 


11.5% 


3.8 MeV 


4° 


50km 


30° 


0.22 


s.5i60 


8.0% 


3.7MeV 


7° 


50km 


40° 


0.20 


s.5i80 


6.1% 


3.6MeV 


8° 


50km 


50° 


0.25 



TABLE II: Properties of the accretion tori at late time. The mass 
of the disk Mdisk (baryon mass outside the excision region), which 
decreases continuously due to accretion onto the BH, is measured 
at imcrgcr + 5ins. Even at late times, all quantities still show 
oscillations of ^ 10%. 



duration of the inspiral varies with the spin of the BH, with 
the component of the BH spin along the orbital angular mo- 
mentum delaying the merger. The merger time tmcigcr, which 
we define as the time at which 50% of the matter has been 
accreted onto the BH, is listed for all cases in Table|T] 

The resulting accretion disk is highly asymmetric, and 
evolves in time. When the disk first forms, around 5ms af- 
ter merger, it creates a torus of matter with its peak surface 
density (baryon density integrated over the height of the disk) 
at r (Einax) ^ 30km and a temperature T ~ 1 — 2McV. Then, 
as matter accretes from the tidal tail and shocks heat the fluid, 
the disk expands quickly. About 10 — 20ms after merger, the 
disk starts to settle into a stable, slowly accreting state. To 
compare the different configurations studied here, we look at 
the properties of this late-time stable configuration, listed in 
Table HI] In Table Hill we give the characteristics of the final 
black hole, as well as the kick velocity, the energy content of 
the emitted gravitational waves, and the peak amplitude of a 
dominant (2,2) mode of the waves. [The (2,2) and (2,-2) are 
the strongest modes, with nearly equal amplitude. See Sec- 
tion lVIBj on the higher modes.] 



Case Mbk/M aBn/MBH i)kick(km/s) Eqw/M rM^l 



sO 


0.97 


0.56 


53 


0.98% 


0.020 


s.5iO 


0.94 


0.77 


60 


0.92% 


0.012 


s.9i0 


0.89 


0.93 


52 


0.95% 


0.009 


s.5i20 


0.95 


0.76 


60 


0.89% 


0.012 


s.5i40 


0.96 


0.74 


61 


0.91% 


0.013 


s.5i60 


0.96 


0.71 


54 


0.95% 


0.014 


s.5i80 


0.97 


0.66 


67 


0.95% 


0.017 



TABLE III: Properties of the post-merger black hole and gravita- 
tional waves. 



A. Effects of spin magnitude 

To test the effects of BH spin magnitude, we compare our 
results for sO, s.5iO, and s.9i0. A comparison of this type 
has already been performed by Etienne et al 11211 . Our cases 
are different, though; unlike them, we do not consider an an- 
tialigned case, but we do push the BH spin to a slightly higher 
level in our run s.9i0. 

Run s.9i0 presented special numerical challenges. In SPEC, 
the singularity and inner horizon of the BH have to be excised 
from the numerical grid while the apparent horizon must re- 
main outside the excision surface. But for nearly extremal 
black holes the region between the inner horizon and the ap- 
parent horizon becomes very narrow. To perform excision in 
such cases, the excision boundary must nearly conform to the 
apparent horizon. We do this by introducing a coordinate map 
in the initial data so that the horizon is initially spherical on the 
pseudospectral grid. We then use our dual frame coordinate- 
control method [34] to fix the location of the horizon through- 
out the whole simulation. For lower spins this is not neces- 
sary, and we only begin to control the horizon location at the 
time of neutron star disruption. That modification excepted, 
case s.9i0 was simulated in exactly the same way as the other 
cases. The deviation between the results at resolutions Res2 
and Res3 is somewhat larger than in the other cases (though 
still quite small). To be more exact, the relative deviation in 
the disk mass between Lev2 and Lev3 is about 9% (~3% of 
the NS mass) for s.9i0 while it was about 5% for sO, and the 
deviation in merger time is about 4% for s.9i0 but only 0.3% 
for sO. The difference indicates that high resolution is needed 
when studying such extreme cases. 

We find that systems with higher aBH/A/sH spiral in more 
slowly: from the same initial separation, sO, s.5iO, and s.9i0 
take roughly 2, 3, and 3.7 orbits, respectively, before NS dis- 
ruption begins. This effect exists in the post-Newtonian treat- 
ment pSy and it has already been seen both in binary black 
hole (e.g. 136]) and BHNS 111211 simulations. Because of the 
prolonged inspiral for the high-spin cases, more angular mo- 
mentum is radiated: O.UAf^ for sO vs 0.141/^ for s.9i0. Ad- 
ditionally, as the BH becomes nearly extremal, increasing the 
spin becomes more and more difficult. This is reflected in the 
final spin of the BH: while aen A^bh increases from to 0.56 
for sO, it only increases from 0.9 to 0.93 for s.9i0. 

We also confirm that the post-merger accretion disk mass 
increases significantly as the magnitude of the aligned BH 
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FIG. 5: Evolution of the total baryonic density outside the hole for 
three different aligned BH spins. 



spin is increased, as shown in Fig. |5] This is in qualita- 
tive agreement with Etienne et al. About 400i\/(llms) after 
merger they find Mb/Mi,,i ^ 20% for an qbyi/Mbyi = 0.75 
system, while for our aBn/A/BH = 0.9 system, we find 
Mb/Mb^i ~35% at a similar time. 



B. Effects of spin orientation 

Most BHNS binaries are expected to have at least a mod- 
erate misalignment between the spin of the BH and the total 
angular momentum lUSll . and this should affect all stages of 
the binary evolution. 

During inspiral, the orbital angular momentum and the BH 
spin precess around the total angular momentum of the sys- 
tem. The evolution of the coordinate components of the BH 
spin for the s.5i80 case is shown in Fig. |6l Over the two 
orbits of inspiral, the spin goes through about a quarter of 
a precession period. The qualitative evolution of the s pin is 
well described by post-Newtonian coiTections (see e.g. 13711 ), 
even though our simulation uses a different gauge choice. 
As for aligned spin, the infall velocity varies between cases: 
the larger the component of the spin aligned with the angu- 
lar momentum, the slower the inspiral. Not too surprisingly, 
we find a monotonic decrease of the merger time with in- 
creasing misalignment angle 0bh, with imerger(0BH, ana) -^ 

tmcrgcr(0,0)for(/)BH->90°. 

The disruption of the star and formation of a disk, shown in 
Fig. |7] proceed somewhat differently from what is observed 
for nonprecessing binaries. As before, the disruption of the 
star is accompanied by the formation of a long tidal tail. But 
because of the inclination of the BH spin, the orbital plane 
of the fluid continues to precess after disruption. Because the 
precession rate varies with the distance to the hole, the tail 
and the disk do not remain in the same plane. While for non- 
precessing binaries matter from the tail falls back within the 
orbital plane of the disk, here matter is added to the disk at 
an angle varying in time. This significantly modifies the na- 
ture of the tail-disk interactions. At small inclination angles 
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FIG. 6: Comparison of the evolution of the Cartesian components 
of the spin for an initial inclination of 80° in our simulation (NR) 
and for the first- and second-order Post-Newtonian expansions (resp. 
IPN and 2PN). The PN values are obtained by integrating the evo- 
lution equations for the spin given in 13711 . using the trajectory and 
current spin of the numerical simulation. 



(<j> ^ 20 ~ 40°), we have a direct collision between the de- 
veloping disk and the tidal tail, while for larger inclinations 
the disk is initially formed of layers of high-density material 
at different angles with respect to the black hole spin. 

The mass of the disk, plotted in Fig. [8] decreases as the in- 
clination of the binary increases. The transition between low 
and high mass disks is continuous, but more rapid at large 
inclinations: for (/)bh < 40°, 10 — 15% of the initial mass 
of the star (^ 0.15 — 0.2AIq) remains either in the tail or 
in the disk 5ms after merger This is roughly similar to the 
disk formed for (/)bh = 0. At higher inclinations, the size of 
the disk drops sharply, to about 5% of the initial mass of the 
star. We expect the disk mass to be even lower for antialigned 
spins (0BH > 90°), though such configurations appear less 
likely to be found in astrophysical systems. These changes in 
the disk mass with the orientation of the BH spin show some 
similarities with the results of Rantsiou et al. |fl6l|, obtained in 
the small mass-ratio limit {q = 1/10) by using a static back- 
ground metric. They found that for a disk to be formed, the 
condition ^bh < 60° has to be satisfied. As our simulations 
use a mass ratio q = 1/3, which is more favorable to the for- 
mation of a disk, it is not too surprising that even high inclina- 
tions leave us with a significant disk; however, the influence 
of inclination remains important for (/)bh > 40°. These fac- 
tors are particularly useful when considering the potential of 
the final remnant to be a short gamma-ray burst progenitor. 
Since the influence of a misaligned BH spin is only felt for 
0BH > 40°, the majority of BHNS binaries can form disks 
about as massive as is predicted by simulations that do not 
take into account the inclination of the BH spin. 

The gravitational wave signal is also significantly affected 
by the value of 0bh- We expand the Newman-Penrose scalar 
^4 extracted at i? = 75A/ using the spin-weighted spheri- 
cal harmonics -2^1™ and choosing the polar axis along the 
initial orbital angular momentum of the binary. The peak am- 
plitude of the dominant (2,2) mode of ^1*4 will increase for 
large inclinations, as could be expected from the results ob- 



tained for aligned spins; here too, a large component of the 
spin along the orbital angular momentum works against large 
^1*4 amplitudes. Additionally, the contribution of subdominant 
modes can become significant at large inclinations. In Fig.|9] 
we show the ratio of the amplitude of the scalar ^4(; m) to 
the amplitude of the (2,2) mode for various {I, ra) modes and 
for 0BH = 20°, 60°. The modes most strongly affected by 
the precession of the binary are the (2,1) and (3,2) modes, 
the first reaching half the amplitude of the dominant mode 
around merger for the s.5i60 simulation. Analytical predic- 
tions for the effect of a precessing trajectory on the modal 
decomposition of the gravitational wave signal have been de- 
rived by Arun et al. 13911 . We find qualitative agreement with 
their results if we assume that the compact objects follow the 
trajectories obtained from our numerical simulations. In par- 
ticular, we note that for precessing binaries, the frequency of 
the (2,1) mode is closer to 2Vt^ot than to Vt^ot- The (2,1) and 
(2,2) modes have similar frequencies, so that the ratio of their 
amplitude computed using the scalar 4*4 is close to the re- 
sult one would obtain by using the gravitational strain h in- 
stead. This will not be true for the (3,3) mode, which has a 
frequency 51(3,3) ~ Sflrot: since ^^ = d^h/dt^, we have 

/i(3,3)//i(2,2) ~ (4/9)*4(3,3)/*4(2,2)- 

In Fig. [TOl we plot the gravitational strain h as observed 
from a distance of 100 Mpc for the simulations s.5i0 and 
s.5i80. The three waveforms correspond to observation points 
whose lines of sight are inclined by 0°, 30° and 60° with re- 
spect to the initial orbital angular momentum of the system. 
Over the short inspiral considered here, the effects of preces- 
sion are relatively small. The main difference visible in these 
waveforms is the slower inspiral experienced by the binary 
with aligned spin. We can also note that in the misaligned 
configuration the star does not disrupt as strongly as in the 
aligned case, causing the cutoff of the wave emission to oc- 
cur at a later time. As a consequence, the cutoff frequency of 
the wave is larger for misaligned spins than for aligned spins. 
This explains why the amplitude of the gravitational strain h 
is comparable for both configurations, while the misaligned 
case showed a significantly larger amplitude when the wave 
was measured using the scalar ^4 (a similar effect occurs if 
both spins are aligned but of different magnitudes). Finally, 
the precession of the orbit causes the gravitational wave emis- 
sion of the misaligned configuration to peak at a non-zero in- 
clination with respect to the initial orbital angular momentum. 
Here, at the time of merger the wave measured at an inclina- 
tion of 30° is slightly larger than at 0°. For these effects to be 
more visible, and in particular for a full precession period to 
be observable, longer simulations are required (^ 10 orbits). 



C. Post-merger accretion disks 

The evolution of the accretion disk over time scales com- 
parable to its expected lifetime is likely to be significantly in- 
fluenced by physical effects that are not taken into account in 
our simulations, mostly the magnetic effects and the impact of 
neutrino cooling. We do not expect the results of our simula- 
tions to accurately represent the details of the late-time evolu- 
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V 





FIG. 7: Evolution of an inclined binary (s.5i80). Top left: Beginning of the simulation, at a separation of 63km. Top right: After 7ms and two 
orbits of inspiral, the star disrupts and most of the matter rapidly accretes onto the black hole. Bottom left: After 1 1ms, the remaining matter 
(~ 8.5% of the initial mass) is split between a developing disk and a tidal tail. Differential precession between the disk and the tail means that 
the disk and the matter falling back from the tail orbit in different planes. Bottom right: After 15.5ms, the disk contains about 4.5% of the 
initial mass of the star. It is still highly inhomogeneous, and slowly expanding. A movie of the whole simulation is available online 13811 . In 
each image, the wired frame shows the "fluid" grid in the z = Q plane (orbital plane at i = 0). 



tion of the disk, but we can nonetheless extract some informa- 
tion regarding the general characteristics of the final remnant. 
To obtain these approximate results, which are summarized in 
Table [III we use the fixed-metric approximation described in 
Sec. lIIB4l starting 5 — 10ms after merger At that time, the 
disk is still expanding, and will typically settle down to a more 
stable quasiequilibrium profile with a relatively low accretion 
rate over about 10ms. 

The coordinate distance between the peak surface density 
of the disk (averaged over all points at a given coordinate ra- 
dius) and the center of the BH shows no strong or monotonic 
dependence on the BH spin. After the initial expansion of the 
disk, variations in the details of the interactions between the 
tidal tail and the accretion disk can lead to different evolutions 



of the density profile. On average, the disks tend to expand 
slightly, while their density decreases because of continued 
accretion onto the black hole. However, neither the tidal tail 
nor the disk are homogeneous, so that the evolution of the den- 
sity profile shows significant oscillations around that average 
behavior. The accretion rate is larger for the more massive 
disks, so that the expected lifetime of the disk is of the same 
order of magnitude for a nonspinning BH (r ^ 75ms) as for 
the highly spinning BH (t ^ 150ms). 

The thermal evolution of the disk does not vary much be- 
tween configurations. The temperature (T)disk rises rapidly 
during the formation of the disk, then stabilizes at about 
3.5MeV in each case — although the average entropy (s)disk 
is about 10% lower for spinning black holes than for sO. For 
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FIG. 8: Evolution of the total baryonic mass outside the hole for 
different initial inclination of the BH spin. 
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FIG. 9: Amplitude of V'(';,m)/V'('2,2) for the 3 modes (2,1), (3,3) and 
(3,2), shown for initial inclinations of the BH spin of 20° and 60° . 
The time f max corresponds to the peak amplitude of the dominant 
(2,2) mode. 



most configurations, the temperature then remains relatively 
stable for the rest of the evolution, with oscillations of order 
10%. The highest spin configuration s.9, however, reaches 
significantly higher temperatures, with T ^ hMeV. Fig.fTTI 
shows the entropy and specific angular momentum profile of 
three of our disks at the end of the simulation. At the final 
time, each of the disks shows an inverted entropy gradient in 
the inner region between the black hole and the radius of max- 
imum surface density. The entropy difference between the in- 
ner edge of the disk and the density maximum is about 10%. 
However, the disk is at least partially stabilized by the strong 
shear in the rotational velocity in these inner regions. 

In Fig. [T2l we show two snapshots of the disk profile for 
the s.5i60 simulation at, respectively, 20ms and 40ms after 
merger. As the disk keeps accreting, the surface density de- 
creases but the profile is otherwise mostly constant. The small 
inverted entropy gradient and the strong positive specific an- 
gular momentum gradient of the inner disk are visible, while 
outside the radius of maximum density the entropy profile is 
mostly constant and the specific angular momentum increases 
more slowly until r k, 100 — 150km. Beyond this radius, mat- 
ter is still in the remnant tidal tail rather than the settled disk. 
The time evolution of these two quantities is extremely small. 
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FIG. 10: Real part of the gravitational strain h for the simula- 
tions S.SiO {top panel) and s.5i80 (lower panel), viewed from differ- 
ent inclinations 6 with respect to the initial orbital angular momen- 
tum. The wave is extracted at r = 75A'f for a mass of the neutron 
star A/ns = 1.4Af0, and assumed to travel as a linear perturba- 
tion up to the observation point located at r = 100 Mpc. Waves 
emitted at the time of merger will reach the radius r = 75Af at 
{t — f merger) ~ 2ms. Over the 2-3 orbits simulated here, the ef- 
fects of the orbital precession — and, in particular, the contribution on 
the second highest mode (2,1), shown in Fig. [9] — remain small. 
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r (km) 



FIG. 11: Upper panel: Specific angular momentum of the accretion 
disk at i = 30ms for runs sO, s.5iO and s.5i60. Note that for the 
smaller spin sO, the disk does not extend farther than r ~ 125km. 
Lower panel: For the same configurations, entropy of the disk aver- 
aged over all points at a given distance from the black hole center. 



the only difference being a smoother profile at late times. The 
disk is relatively thick, with H/r ~ 0.2 at all radii within the 
disk, a value that remains constant from a few milliseconds 
after merger to the end of the simulation. 

For inclined disks, we also measure the tilt /? and twist 7 of 
the disk, as defined by Eq. ( l22b . For large inclinations (s.5i60 
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FIG. 12: Profile of the disk forming in the simulation s.5i60 at t 
20ins (upper panel) and t — 40ms (lower panel). 



and s.5i80), the tilt angle between the disk angular momen- 
tum and the orientation of the BH spin varies slowly during 
the evolution of the disk, at least in the higher density regions. 
Fig. [13] shows the tilt and twist profiles 20, 30 and 40ms af- 
ter merger for simulation s.5i60. In the inner part of the disk 
(r < 50km), the tilt decreases from 15° to 7° over those 20ms. 
The precession of the disk is more significant: we observe 
a variation of the twist of about 100° over the same period. 
However, the precession rate is not constant at all radii and 
changes in time; the disk has not reached a state in which it 
precesses at a constant rate as one solid body. For smaller 
inclinations, relative variations in the tilt are larger. The in- 
clination of the disk decreases at late-time to /3 < 5°, and no 
global precession is observed. 

An inclined disk showing some similarities with the re- 
sults of our most inclined simulations (s.5i60 and s.5i80) 
was evolved in the presence of magnetic fields by Fragile et 
al. ll40ll . Even though the two simulations vary greatly in their 
initial conditions — Fragile et al. start their simulation from a 
torus of matter with peak density at r = 25A/ {^ 200km) — 
the thickness and inclination of the disks are equivalent and 
some results from ll40ll could apply to the late time behavior 
of our disks. In Fragile et al., the inner disk is warped by 
the gravitomagnetic torque of the black hole, leading to larger 
tilts at lower radii. The same torque leads to a precession of 
the disk over a period of about 4s (for a black hole remnant 
of final mass A/bh ~ 5.6A/o). This last effect is, however, 
acting over time scales longer than the lifetime of the disk 
formed in BHNS mergers, so that it seems unlikely that our 
disk would have time to reach the steady precession described 
in [40l:]. The magneto-rotational instability (MRI), on the other 
hand, appears to develop over roughly one orbital time scale, 
or about 20ms for the initial configuration chosen in 114011 and 
a central black hole of mass A/bh = 5.6A/0. The rise of the 
MRI might be even faster for a disk more similar to the results 
of our simulations, as the peak of the density profile is signifi- 
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FIG. 13: Upper panel: Tilt profile of the accretion disk formed 
in simulation s.5i60, as obtained from evolutions on a static back- 
ground metric. The inclination of the disk decreases in time, with 
/3 ~ 10 — 15° at the time of disk formation, but /3 ~ 5 — 7° towards 
the end of the simulation. 

Lower panel: Twist profile for the same configuration. Over the 
20ms of evolution, the disk goes through more than one fourth of 
a precession period. 



cantly closer to the black hole in our disks than in 114011 . and the 
evolution time scale is thus shorter: the orbital period of circu- 
lar orbits at r = 50km is about 5ms. The MRI should have a 
strong influence on the redistribution of angular momentum in 
the fluid, and therefore on the accretion rate. However, the ac- 
cretion rate is also influenced by the presence of shocks in the 
disk (see e.g. Fragile and Blaes |41] for shocks in tilted disks 
similar to ll40ll ). This means that interactions between the disk 
and the matter falling back in the tidal tail are also likely to 
play an important role in the determination of the lifetime of 
the disk. Thus, both the magnetic effects and realistic initial 
conditions are required to accurately predict the lifetime of the 
disks resulting from BHNS mergers. 



VII. CONCLUSIONS 

Astrophysical BHNS binaries are expected to have BH 
spins that are not aligned with the orbital angular momen- 
tum of the system. We performed here the first fully gen- 
eral relativistic simulations of BHNS systems with precess- 
ing orbits. We find that for realistic inclinations of the BH 
spin with respect to the initial orbital angular momentum 
(0BH = — 80°), a mass ratio of 1:3, and a moderate black 
hole spin aBn/^lBH = 0.5, the mass of the disk varies by 
about a factor of 2. More important, the inclination of the 
spin seems to have a significant impact on the final remnant 
only for (J>bk > 40°. According to population synthesis mod- 
els by Belczynski et al. jTsil . this means that, for binaries with 
initial spin of 0.5, half of the systems would have disks nearly 
as massive as if the BH spin was aligned with the orbital angu- 
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lar momentum of the system. This confirms the relevance of 
the aligned BHNS studies previously undertaken by ourselves 
and by other groups. At late times, the inclination of the disks 
formed in these precessing systems remain relatively modest 
(/3 < 15°). Most of the angular momentum of the system is in 
the orbital motion of the binary, which precesses around the 
total angular momentum of the system at a small misalign- 
ment angle (- 5 - 20° for 0bh == 20 - 80°). The black 
hole spin axis itself is inclined at a larger angle to the total an- 
gular momentum, but its misalignment decreases as the hole 
accretes matter from the disrupted neutron star This suggests 
that more inclined disks could be observed for larger black 
hole spins or more extreme mass ratios. Here, our most in- 
clined binaries have an average tilt [3 ^ 10°. From the results 
of Fragile et al. [4(11, we would expect those disks to precess 
as one solid body around the black hole — but only over time 
scales far longer than the expected lifetime of our post-merger 
disks (Tpi-cc ~ 4s >> Tacc ~ lOOms). 

For spins aligned with the orbital angular momentum, our 
study shows qualitative agreement with previous results by 
Etienne et al. II12I1 . As expected, large spins favor the for- 
mation of a massive disk. By studying a higher initial spin 
(aBn/^-Z^BH = 0.9), we also show that large disks of mass 
Mdisk ~ O.5A/0 can be obtained for A/bh == 4:.2Mq. 

All the cases studied here produce post-merger systems that 
are promising SGRB central engines. Despite large differ- 
ences in the disk masses, the lifetime of the system seems 
mostly independent of the black hole spin; we find an accre- 
tion time scale Tacc ~ 75 — 150ms for all cases. The disks 
have a peak surface density located about 50km away from 
the hole, and extend about twice as far They are always thick 
{H/r w 0.2), hot ((T) = 3 - 5MeV), and accreting at a super- 
Eddington rate (Af = 0.5 — 5Mq/s). All simulations also 
show the presence of a baryon-free region, at least at densities 
above the threshold at which atmospheric corrections begin 
to have an impact (~' lO^g/cm'^). This region covers a cone 



with an opening angle of 30 — 50° around the axis of the black 
hole spin, except for the high-spin configuration for which the 
disk is significantly closer to the BH, and the opening angle 
varies within the range 5 — 10°. Such a region is required if 
relativistic jets are to be launched. 

Accretion continues throughout the disk evolution, but the 
thermal and rotational profiles do seem to stabilize. The spe- 
cific angular momentum of the disks increases with radius, so 
they are not subject to the Rayleigh instability. The angular 
velocity, however, decreases with radius, so these disks are 
subject to the magneto-rotational instability (MRl), an effect 
not included in our simulations. 

The late-time behavior of the black hole-accretion disk sys- 
tem is critical if we want to understand the potential of BHNS 
mergers as progenitors for short gamma-ray bursts. Currently, 
the measurement of the properties of the disk and their evo- 
lution in time suffers from the limitations of our simulations. 
The general characteristics of the disk can be obtained, but a 
more detailed evolution would certainly require the inclusion 
of magnetic fields and neutrino radiation. These effects will 
be added to our evolutions in the near future. 
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